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ABSTRACT 

We calculate the secondary anisotropics in the CMB produced by inhomogeneous 
reionization from simulations in which the effects of radiative and stellar feedback 
effects on galaxy formation have been included. This allows to self-consistently 
determine the beginning (zj ~ 30), the duration (5z ~ 20) and the (nonlinear) 
evolution of the reionization process for a critical density CDM model. In addition, 
from the simulated spatial distribution of ionized regions, we are able to calculate the 
evolution of the two-point ionization correlation function, C^, and obtain the power 
spectrum of the anisotropics, Ci, in the range 5000 < i < 10^. The power spectrum 
has a broad maximum around £ ~ 30000, where it reaches the value 2 x 10~^^. We 
also show that the angular correlation function C{9) is not Gaussian, but at separation 
ang les e < 10"^ rad it can be approximated by a modified Lorentzian shape; at larger 
separations an anticorrelation signal is predicted. Detection of signals as above will be 
possible with future mm-wavelength interferometers like ALMA, which appears as an 
optimum instrument to search for signatures of inhomogeneous reionization. 



1. INTRODUCTION 

At z ~ 1100 the intergalactic medium recombined and remained neutral until the first 
sources of ionizing radiation form and begin to reionize it. Current models of cosmic structure 
formation predict that the first collapsed, luminous objects should have formed at redshift 
z ~ 30. This conclusion is reached by requiring that the cooling time tc of the gas be shorter 
than the Hubble time tn at the formation epoch. The ionizing fiux from these objects creates 
cosmological HII regions in the surrounding IGM, whose sizes are much smaller than their typical 
interdistance (Ciardi, Ferrara & Abel 1999). This implies that reionization passed through a 
highly inhomogeneous phase which ended only when the individual HII regions overlapped, or 
stated differently, when reionization was complete. There is essentially no direct observational 
test of both the beginning and the duration of the reionization process. The most stringent 
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upper limits on the reionization redshift and the Thomson optical depth of the cosmic medium, 
derived from data on "linear" anisotropics at £ < 10^, are quite model-dependent (De Bernardis et 
al. 1997, Griffiths et al. 1999). 

However, if reionization was indeed inhomogeneous, it should have left an inprint in the CMB. 
In the homogeneous case, small-scale secondary anisotropies generated by Dopplcr effect from 
Thomson scattering off IGM electrons would be erased due to potential flow cancellation effects. 
In the inhomogeneous case, the modulation of the ionization fraction, playing a similar role as 
the density modulation for the nonlinear Vishniac effect, prevents such cancellation leading to 
anisotropies at sub-degree scales. A few papers have recently tackled the calculation of secondary 
anisotropics produced by inhomogeneous reionization (Aghanim et al. 1996 [see also Erratum 
(1999)], Knox, Scoccimarro & Dodelson 1998, Gruzinov & Hu 1998, Peebles &: Juszkiewicz 1998, 
Haiman &: Knox 1999, Hu 1999). These models are based on simple assumptions concerning the 
evolution of the mean ionization level, the typical size of the ionized patches and, in some cases, 
their spatial and redshift correlation. 

Our aim here is to improve the modelling by self-consistently calculating secondary 
anisotropies from dedicated reionization simulations presented in Ciardi, Ferrara, Governato &; 
Jenkins (1999, CFGJ), in which the effects of radiative and stellar feedback effects on galaxy 
formation have been included. Similar reionization simulations have been performed also by 
Gnedin (1998, 1999). This allows us to calculate the history and correlation properties of 
reionization to an unprecedented detail level. 

2. INHOMOGENEOUS REIONIZATION 

The predictions for the CMB anisotropies presented in this paper are based on the study of 
inhomogeneous reionization (IHR) presented by CFGJ. Although the detailed description of the 
model can be found in that paper, it is useful to recall here its essential features and the results 
relevant for the present work. 

The reionization process is studied in a critical density CDM universe (JIq = 1, h = 0.5 
with as = 0.6 at 2; = and Uj, = 0.06) as due to stellar sources, including Population III 
objects. The spatial distribution of the sources is obtained from high-resolution numerical N-body 
simulations within a periodic box of comoving length L = 2.55h~^ Mpc. The source properties 
are calculated taking into account a self-consistent treatment of both radiative (i.e. ionizing and 
H2-photodissociating photons) and stellar {i.e. SN explosions) feedbacks regulated by massive 
stars. This allows, in particular, to derive the spatial distribution of the ionized regions at various 
cosmic epochs and the evolution of the mean H ionization. In brief, there are two main free 
parameters in the simulations: (i) the fraction of total baryons converted into stars and 
(ii) the escape fraction of ionizing photons, fesa from a given galaxy; a critical discussion of 
these parameters is given in CFGJ. Four different combinations of these parameters in runs A 
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(/^ = 0.012, /esc = 0.2), B (0.004, 0.2), C (0.15, 0.2) and D (0.012, 0.1) have been explored. Run 

A gives the best agreement between the derived evolution of the cosmic star formation rate and 
the experimentally deduced one at z ^ 4 (Steidel et al. 1998, CFG J). Therefore, here we only 
present results for this case. 

The topological structure of ionization at z = 19.8 for run A is shown in Fig. 1. The proper 
size of ionized bubbles is approximately in the range 1-20 kpc at this epoch, as only relatively 
small, and hence faint, objects have collapsed; they grow in number and volume with cosmic time.^ 
A more global view of the reionization process is given in Fig. 2, where the redshift evolution of 
the volume-averaged mean ionization fraction, x(z) is shown for the four different runs. Except 
for run C (high star formation efficiency), when reionization is complete at z ~ 15, primordial 
galaxies are able to reionize the IGM at a redshift z ~ 10. Note that reionization begins at Zi ~ 30 
when the conditions of the cosmic medium allow the first Pop III objects to start to form, and 
it is then completed after a redshift interval Sz 20 thanks to the contribution of larger objects 
with mass ^ lO^M©. 



3. SECONDARY ANISOTROPIES 

To calculate the secondary anisotropics produced by IHR, we use the same method outlined 
by Knox et al. (1998) and Gruzinov Sz Hu (1998). The solution of the Boltzmann equation for the 
present value of the perturbation of the photon temperature A = ST/T can be written as 

Ao(7) = To / -jx{^, V)l ■ v(x), (1) 

where Hq = lOO/i km s~^ Mpc~^ is the Hubble constant, 77 = 2/Hq{1 + z)^/^ is the conformal time. 
To = nefiCFTVoc = 0.137ri(,/t//Lt with ne,o being the present free electron density, ut the Thomson 

cross section, and ^ = 0.59 the mean molecular weight of a cosmological mixture of ionized H 
and He. The quantity x(X) 'H) the ionization fraction calculated at position x = 7(??o ~ v) ''■i^d 
conformal time 77; v(x) is the peculiar velocity today in units of c. The two-point correlation 
function due to ionization is then 

C{e) = (Ao(7l)Ao(72))|cose=7l-72, (2) 

or, using the expression for Ao(7) above, 

c{e) = Tl /'^ ^c,{m,m,o)c^{m,m,e), (3) 

Jo 111 •'0 V2 

where Cy{r]i,r]2,9) = (71 • v(xi)72 • v(x2)) and C^(r]i,r]2,9) = (x(xi, r/i)x(x2, 7/2)) (we have 
implicitly assumed that the velocity and ionization are independent fields) . The velocity correlation 
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function depends on the matter power spectrum P{k). Groth et al. (1989) have shown that Cy 
can be written as 

C,{Vi,V2,0) = PuHir) + Si2^{r), (4) 

where Pi2 = (71 • r)(72 • r)/r, 6*12 = 7i • 72 ^ P12 and r = xi — X2. The functions 11 (r) and E(r) 
can be written (for a vanishing cosmological constant) as 

n(r) = iHin'-''/2Tr'') / dkP{k)jo{kr) - 2S(r), (5) 
Jo 

/•CO 

S(r) = {H^n^-'^ /2Tr'^) / dkP{k){kr)-^ ji{kr). (6) 
We have adopted a CDM power spectrum of fluctuations as in Efstathiou, Bond & White (1992): 

P{k) = Ak[l+[Bk + {Ckfl^ + {DkfYY'^^'' ' (7) 

where u = 1.13, B = 6.4 /i"^ Mpc, C = 3.0 /i"^ Mpc and D = 1.7 ft,"^ ^pc. The power spectrum 
and its normalization are consistent with the above numerical simulations. 

The most delicate point of the calculation of reionization anisotropies is probably the form 
of C^^. Previous studies (see the Introduction) have often computed this function assuming a 
superposition of ionized patches of equal size R. In addition, the initial redshift of reionization Zj, 
its duration 5z, and the evolution of the mean ionization fraction xiv) to be postulated there. 
Finally, the spatial correlation of the ionizing sources is also crucial (for a discussion see Oh 1998) 
as it might induce analogous correlations in the ionization field. In their analytic treatment Knox 
et al. 1998 have been forced to postulate an ad-hoc functional form for this quantity. Thanks 
to the IHR simulations presented above, we are instead able to calculate without further 
assumptions as explained in the next Section. 



4. THE IONIZATION CORRELATION FUNCTION 

To calculate we proceed as follows. Our simulations provide us with 11 data 
boxes containing the spatial information on the ionization field for the following redshifts 
z = 8.3, 10.9, 11.4, 14.3, 15.4, 16.5, 18.0, 19.8, 22.1, 25.3, and 29.6. As an example, we have shown in 
Fig. 1 the one corresponding to z = 19.8. We label each epoch with an index n = 1, 9, where 
n = 1 (n = 9) corresponds to z = 10.9 (z = 25.3). For each box we select a slice through its 
center, thus obtaining a N x N = 256^ array with the values of the ionization fraction at each 
cell point. The first and the last slice have X = 1 (fully ionized) and x = (completely neutral), 
respectively. Therefore they do not contribute to anisotropies from inhomogeneous reionization 
and can be disregarded in the calculation. For each couple of slices with redshifts Zm we then 
calculate the total two-point correlation function as 

cr(^,.n,.„.)=Ei:^^^^^%^. (8) 
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The {1/2)N{N — 1) components of C^"* (for each couple n,m) thus obtained are binned according 
to the separation angle 9 between the two lines of sight passing through the cells i and j. Hence, 
9 runs from 9jnin = 2.0 x 10~^ rad, the angle subtended by the cell proper size at the highest 
redshift of the simulation, to 9max = 8.4 x 10"'^ rad, corresponding to proper size V2L at the 
lowest redshift considered; the angular distance used is the standard Priedmann one. Ng is then 
the number of values falling in [9,9 + 9min]- We repeat the procedure for a random slice realization 
with the same volume-averaged mean of the ionization, thus obtaining Cr{9). Assuming Poisson 
noise, the associated statistical error on the correlation functions is C{9)/\/N^. 

The correlation function due to ionized patches {i.e. the inhomogeneous part) can then be 
written as: 

'^^^ ^"*) (^1 ^rri)- (9) 

This function correctly describes reionization effects due to ionized patches: for example, at the 
lowest redshift of our simulations, where the reionization is complete, Cy.{9) = 0. 

In Fig. 3 we show the autocorrelation functions C"" along with the corresponding errors; the 



largest error in the curves is 1.9% The behavior of C^^ is rather flat for small angles (0 ^ 10 ® rad) 



X 

^x 

and decreases rapidly for larger separations; for even larger values of 9, becomes negative, 
implying that the ionization regions are anticorrelated. The position of the first zero, 9o, 
provides an estimate of the typical size of the ionized regions (including overlapping) at the 
corresponding redshift through the relation R ~ 0.23(^o/10~^rad)/i~-^ Mpc. The amplitude of the 
anticorrelation signal is roughly proportional to the correlation one. For small separations C^*" 
can be approximated by a modified Lorentzian function of the type 
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1 + [(^ - emm)/r]" ' 

where A^ax is the maximum amplitude at 9 = 9min, and F is the width, a is the power index. 
The values of these fitting parameters for the autocorrelation functions are given in Tab. 1. 

The anticorrelated regions cannot be described in terms of such a simple assumption. 
Therefore the shape is different from the Gaussian one proposed by Gruzinov & Hu 1998) with 
redshift-dependent amplitude and standard deviation. The autocorrelation functions for any n are 
always larger than those obtained by correlating slices at different redshifts (m ^ n), although the 
shape remains similar. We have also checked that the results do not depend on the particular slice 
used provided it is located not too far from the center of the simulation box. Off-center planes 
might lead to sensibly different amplitudes of (« ±20%) due to border effects, only if their 
position is closer than 0.04L to the box walls. 



5. RESULTS AND CONCLUSIONS 



By using the derived ionization correlation function, we can calculate the angular power 
spectrum produced by the inhomogeneous part of the reionization from eq. 2. In Fig. 4 we show 
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the behavior of the two-point correlation function, C{6), and of the angular power spectrum Q 
which, for ^ S> 1, can be approximated by 

roo 

Q = 27r / deejQ{ie)c{e). 

Jo 

The correlation function smoothly decreases from its maximum C(0) 5 x 10^^^ and becomes 
negative at 6 = 1.7 x 10~^ rad as a result of the anticorrelation discussed above. Beyond the 
second null point ~ 5 x 10^^ rad) it rapidly becomes more noisy due the decreasing number of 
data points in the simulation at these large separations. In addition to the statistical error, this 
region is also affected by systematic errors. These are mainly due to the finite size of the box 
which does not allow to take into account ionized regions produced by sources located just outside 
the simulation volume. 

To minimize the of this effect we have cut the spectrum beyond separations that are larger 
than the minimum difference between the angular size of the box and the typical angular size 
of the ionized regions (see previous Section). This procedure should yield a C{9) reliably 
unaffected by systematic errors below ~ 6 x 10~^ rad: we therefore limit our analysis to the 
angular interval Oxnin 

The reported angular power spectrum therefore extends in the range 5000 ^ £ ^ 10^ because 
of the relation £ = tt/9. The spectrum shows an absolute maximum around i ~ 30000, where it 
has the value 2 x 10~^^. Beyond £ ~ 10^ ringings appear in Fig. 4 which are caused by the abrupt 
cutoff of C{9) at 9c, but up to ^ 10^ the behaviour of angular spectrum is not very far from 
flatness. We have calculated a running average Se of the spectrum (also shown in Fig. 4) through 
the formula 

3e/2 

= ^ E ^^A' (11) 

x=e/2 

This allows to get rid of the ringings (which have essentially zero average) and to appreciate the 
decrease to zero of the spectum above £ = 2 x 10"''. Our results match very closely those obtained 
by Knox et al. (1998). In particular, the best agreement is found with their correlated-source 
model corresponding to Zi = 35, Sz = 4. On the other hand, our numerically derived values 
are Zi 30, Sz Ri 15 (also notice that our evolution of xiz) is nonlinear). Therefore the close 
agreement can be regarded as rather accidental, given also the various assumptions made in that 
work; we recall that we have derived here the spatial distribution and evolution of the ionization 
fraction from numerical simulations in which several feedback effects regulating the formation of 
the ionizing sources have been included. Also shown for comparison in Fig. 4 are the primary 
angular spectra corresponding to a CDM model with the same parameters adopted here, in which 
the effects of a homogeneous reionization with total r = 0.177 (as derived from our simulations) 
are either neglected or includuded. These spectra have been obtained by running CMBFAST 
(Seljak & Zaldarriaga 1996). The maximum value of our angular spectrum is comparable to, but 
higher than the most recent estimates of the Vishniac signal (Jaffe &; Kamionkowski 1998, Hu 
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1999). For a cluster-normalized CDM model (which is thereby directly comparable with ours), 
Jaffe & Kamionkowski (1998) find a peak signal ~ 8 x 10~^^, which is however located at lower 
multipole numbers, around H ~ 5000. We also note that our predicted maximum value is typically 
smaller than the peak anisotropy in un uncorrelated-source models with large ionized bubbles. 
For instance it is about 3 times smaller than the estimate given by Gruzinov &; Hu (1998), for 
their favorite model with Zi = 30, (^2; = 5, i? = 3 Mpc [the difference is smaller if we consider more 
refined computation of Knox et al (1999)]. For such uncorrelated-source models however a large 
bubble size implies that the peak is confined to ^ ^ 5000. In any case, all of these models predict 
that the signal from IHR should dominate both the primary and the Vishniac one for I ^ 5000. 

The detection of this signal would be of great importance for the study of early galaxy 
formation, the IGM and the nature of the reionization sources. If reionization is predominantly 
produced by stellar type objects, as assumed here, the relatively small, sub-Mpc sizes of the 
ionized patches and their correlation properties, tend to shift the peak of the power spectrum 
towards larger multipole values, £ 30000. This detection will require the use of the next 
generation of millimiter wavelength interferometers like ALMA ^.This instrument is expected to 
reach sensitivities of 2 fiK in one hour and reach I ~ 10^, thus appearing as a perfect instrument 
to search for signatures of IHR, and, in general, for the above type of studies. Coupled with 
theoretical predictions for different cosmological models and simulations exploring a wide range of 
variation of the key parameters, as /b* and fesc, these experiments can give a comprehensive and 
clean view of how reionization proceeded in the universe. 

Instruments nearer in the future as MAP and PLANCK could at least clarify and constrain 
the role, if any, played by quasars in the reionization process. Mini-quasars have been advocated 
by some authors (Haiman, Madau Sz Loeb 1999) as possible sources for reionization; in this 
scenario, we would expect larger patch sizes and a spectrum peaked at lower values of i, as in 
the uncorrelated-source models with R of order a few Mpc. However, uncertainties might remain 
related to the duty-cycle of such objects, if they exist. 

IHR might in principle affect the determination of cosmological parameters planned with 
future experiments. Our calculations do not allow us to quantify the uncertainty introduced 
because the finite size of the simulation box prevents extension of our results below £ 5000, 
where the primary spectrum should be already more than one order of magnitude smaller as 
shown by Fig. 4. However, since from our results the angular spectrum appears to be already 
declining for decreasing I for I < 10^, a comparison with the results of Knox et al. (1998) seems 
to imply that the impact of IHR on PLANCK may somewhat less than in their analytical model. 
Those authors also conclude that the very small scales affected by IHR do not sensibly alter the 
parameter determination obtained by MAP; our results further support their conclusion. 

Although our calculations do not include the polarization angular spectrum Cpi, a rough 
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estimate can be obtained extrapolating the analytic results of Hu (1999) to the present case. 
According to this author CpijCn ~ 10~^ (Qrms/^^rms)^) where the ratio of rms quadrupole to rms 
velocity field is evaluated at z = assuming bzjzi ^ 1. Although the simulations show that 
bzjzi K, 1, we can nevertheless infer that CpijCi k. 10~^. Since for r = 0.177 the "linear" angular 
spectrum has CpijCi ^ 10~^ for large i, we conclude that the impact of IHR on cosmological 
parameter determination by PLANCK will be even less when polarization data are taken into 
account. The joint exploitation of temperature and polarization power spectra (Efstathiou & 
Bond 1998, Kinney 1998) is in fact required for the removal of some cosmological parameter 
degeneracies. 

We would like to thank R. Scoccimarro for useful suggestions and P. de Bernardis, S. Masi, 
and V. Natale for stimulating discussions. This work is partly supported by a MURST-COFIN98 
(AF) and Agenzia Spaziale Italiana-ASI (RF) grants. 
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1.7 X 10-5 


1.7 



Table 1 : Autocorrelation function fitting parameters 



Fig. 1 



Figure available at http://www.arcetri.astro.it/~ferrara/reion.html 

, — Snapshot of the distribution of ionized regions in the simulation box at redshift z = 19.8 for run A 
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Fig. 2. — Evolution of the volume- averaged mean ionization, x(.^)y for the four runs A, B, C, D whose parameters 
are given in the text. 
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Fig. 3. — Evolution of the ionization autocorrelation function C^"{9) with n — 1, .., 9, where n labels the redshift 
as explained in the text. Also shown are the statistical errors (shaded areas). 
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Fig. 4. — Top: Angular correlation function C{9); Bottom: Power spectrum d of secondary anisotropies produced 
by IHR; superposed is the running spectrum average Se- The curves on the left are the primary spectra in which 
homogenous reionization effects are neglected (solid) or included (dashed) with a reionization optical depth r — 0.177. 



